home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / mathpack / sgefa.f < prev    next >
Text File  |  1989-08-17  |  3KB  |  104 lines
  1.       subroutine sgefa(a,lda,n,ipvt,info)                               
  2.       integer lda,n,ipvt(1),info
  3.       real a(lda,1)
  4. c
  5. c     sgefa factors a real matrix by gaussian elimination.
  6. c
  7. c     sgefa is usually called by sgeco, but it can be called
  8. c     directly with a saving in time if  rcond  is not needed.
  9. c     (time for sgeco) = (1 + 9/n)*(time for sgefa) .
  10. c
  11. c     on entry
  12. c
  13. c        a       real(lda, n)
  14. c                the matrix to be factored.
  15. c
  16. c        lda     integer
  17. c                the leading dimension of the array  a .
  18. c
  19. c        n       integer
  20. c                the order of the matrix  a .
  21. c
  22. c     on return
  23. c
  24. c        a       an upper triangular matrix and the multipliers
  25. c                which were used to obtain it.
  26. c                the factorization can be written  a = l*u  where
  27. c                l  is a product of permutation and unit lower
  28. c                triangular matrices and  u  is upper triangular.
  29. c
  30. c        ipvt    integer(n)
  31. c                an integer vector of pivot indices.
  32. c
  33. c        info    integer
  34. c                = 0  normal value.
  35. c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
  36. c                     condition for this subroutine, but it does
  37. c                     indicate that sgesl or sgedi will divide by zero
  38. c                     if called.  use  rcond  in sgeco for a reliable
  39. c                     indication of singularity.
  40. c
  41. c     linpack. this version dated 08/14/78 .
  42. c     cleve moler, university of new mexico, argonne national lab.
  43. c
  44. c     subroutines and functions
  45. c
  46. c     blas saxpy,sscal,isamax
  47. c
  48. c     internal variables
  49. c
  50.       real t
  51.       integer isamax,j,k,kp1,l,nm1
  52. c
  53. c
  54. c     gaussian elimination with partial pivoting
  55. c
  56.       info = 0
  57.       nm1 = n - 1
  58.       if (nm1 .lt. 1) go to 70
  59.       do 60 k = 1, nm1
  60.          kp1 = k + 1
  61. c
  62. c        find l = pivot index
  63. c
  64.          l = isamax(n-k+1,a(k,k),1) + k - 1
  65.          ipvt(k) = l
  66. c
  67. c        zero pivot implies this column already triangularized
  68. c
  69.          if (a(l,k) .eq. 0.0e0) go to 40
  70. c
  71. c           interchange if necessary
  72. c
  73.             if (l .eq. k) go to 10
  74.                t = a(l,k)
  75.                a(l,k) = a(k,k)
  76.                a(k,k) = t
  77.    10       continue
  78. c
  79. c           compute multipliers
  80. c
  81.             t = -1.0e0/a(k,k)
  82.             call sscal(n-k,t,a(k+1,k),1)
  83. c
  84. c           row elimination with column indexing
  85. c
  86.             do 30 j = kp1, n
  87.                t = a(l,j)
  88.                if (l .eq. k) go to 20
  89.                   a(l,j) = a(k,j)
  90.                   a(k,j) = t
  91.    20          continue
  92.                call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
  93.    30       continue
  94.          go to 50
  95.    40    continue
  96.             info = k
  97.    50    continue
  98.    60 continue
  99.    70 continue
  100.       ipvt(n) = n
  101.       if (a(n,n) .eq. 0.0e0) info = n
  102.       return
  103.       end
  104.